rm(list=ls())
library(readstata13)
library(car)
library(ggplot2)
library(nnet)
library(MASS)
library(nnetpredint)
library(stargazer)
library(ordinal)
library(glm.predict)
library(plyr)
library(stringr)
library(zoo)
setwd("C:/Users/gschuma1/Dropbox/Papers/Opinion_polls/")
exp.data <- read.dta13("Opinion_Schumacher.dta", convert.factors=FALSE)


# Do some recoding --------------------------------------------------------
party <- recode(exp.data$party_combined,"1='Left Party';
                                         2='Social Democrats';
                                         3='Center Party';
                                         4='Liberal Party';
                                         5='Conservatives';
                                         6='Christian Democrats';
                                         7='Green Party'; 
                                         8='Sweden Democrats';
                                         9='Feminist Initiative';
                                        10='Other'")

exp.data$party_close_rec <- recode(exp.data$party_close_h14,"1=0;
                                                             2=1;
                                                             3=1;
                                                             4=1")

exp.data$change_econ_rec <- recode(exp.data$Change_economy, "3=0;
                                                             2=1;
                                                             4=1;
                                                             5=2;
                                                             1=2")

exp.data$change_welfare_rec <- recode(exp.data$Change_welfare, "3=0;
                                                                2=1;
                                                                4=1;
                                                                5=2;
                                                                1=2")
exp.data$Change_policy_overall <- exp.data$Change_policy_overall-1
exp.data$policy_change <- recode(exp.data$Change_policy_overall,"0=2;
                                                                 2=0")
exp.data$policy_change2 <- recode(exp.data$policy_change,"2=1")
exp.data$econ_rec2 <- recode(exp.data$change_econ_rec,"2=1")
exp.data$welfare_rec2 <- recode(exp.data$change_welfare_rec,"2=1")

exp.data$change_welfare2 <- recode(exp.data$Change_welfare, "5=2;
                                                             4=2;
                                                             3=1;
                                                             2=0;
                                                             1=0")


exp.data$change_economy2 <- recode(exp.data$Change_economy, "5=2;
                                                             4=2;
                                                             3=1;
                                                             2=0;
                                                             1=0")
exp.data$Change_strategy2 <- recode(exp.data$Change_strategy,"1=0;
                                                              2=0;
                                                              3=1;
                                                              4=2;
                                                              5=2")

exp.data$gov <- ifelse(party=="Social Democrats",1,
                ifelse(party=="Green Party",1,0))
exp.data$seem.to.lose <- ifelse(party%in%c("Social Democrats", "Green Party", "Liberal Party",
                                           "Christian Democrats", "Feminist Initiative"),1,0)
exp.data$party <- party

# Was the reference category treatment long ago or not
exp.data$n.months <- ifelse(exp.data$party_combined==1 & exp.data$Treatment_polls==0, 25,
                     ifelse(exp.data$party_combined==2 & exp.data$Treatment_polls==0, 27,       
                     ifelse(exp.data$party_combined==3 & exp.data$Treatment_polls==0, 8,
                     ifelse(exp.data$party_combined==4 & exp.data$Treatment_polls==0, 25,
                     ifelse(exp.data$party_combined==5 & exp.data$Treatment_polls==0, 37,
                     ifelse(exp.data$party_combined==6 & exp.data$Treatment_polls==0, 16,
                     ifelse(exp.data$party_combined==7 & exp.data$Treatment_polls==0, 7,
                     ifelse(exp.data$party_combined==8 & exp.data$Treatment_polls==0, 6,
                     ifelse(exp.data$party_combined==1 & exp.data$Treatment_polls==1, 18,
                     ifelse(exp.data$party_combined==2 & exp.data$Treatment_polls==1, 16,       
                     ifelse(exp.data$party_combined==3 & exp.data$Treatment_polls==1, 25,
                     ifelse(exp.data$party_combined==4 & exp.data$Treatment_polls==1, 16,
                     ifelse(exp.data$party_combined==5 & exp.data$Treatment_polls==1, 26,
                     ifelse(exp.data$party_combined==7 & exp.data$Treatment_polls==1, 1,
                     ifelse(exp.data$party_combined==8 & exp.data$Treatment_polls==1, 15,NA   
                            )))))))))))))))

exp.data$gain.dif <- ifelse(exp.data$party_combined==1 & exp.data$Treatment_polls==0, (6.8 - 8) / 6.8,
                     ifelse(exp.data$party_combined==2 & exp.data$Treatment_polls==0, (29.5-34.3) / 29.5,       
                     ifelse(exp.data$party_combined==3 & exp.data$Treatment_polls==0, (6.1-7.4) / 6.1,
                     ifelse(exp.data$party_combined==4 & exp.data$Treatment_polls==0, (5.4-6.3) / 5.4,
                     ifelse(exp.data$party_combined==5 & exp.data$Treatment_polls==0, (24.7-28.3) / 24.7,
                     ifelse(exp.data$party_combined==6 & exp.data$Treatment_polls==0, (3.1 - 4.3) / 3.1,
                     ifelse(exp.data$party_combined==7 & exp.data$Treatment_polls==0, (4.7-6.3) / 4.7,
                     ifelse(exp.data$party_combined==8 & exp.data$Treatment_polls==0, (17.3 - 20) / 17.3,
                     ifelse(exp.data$party_combined==1 & exp.data$Treatment_polls==1, (6.8 - 5.6) / 6.8,
                     ifelse(exp.data$party_combined==2 & exp.data$Treatment_polls==1, (29.5-24.8) / 29.5,       
                     ifelse(exp.data$party_combined==3 & exp.data$Treatment_polls==1, (6.1-4.6) / 6.1,
                     ifelse(exp.data$party_combined==4 & exp.data$Treatment_polls==1, (5.4-4.8) / 5.4,
                     ifelse(exp.data$party_combined==5 & exp.data$Treatment_polls==1, (24.7-20.9) / 24.7,
                     ifelse(exp.data$party_combined==7 & exp.data$Treatment_polls==1, (4.7-3.5) / 4.7,
                     ifelse(exp.data$party_combined==8 & exp.data$Treatment_polls==1, (17.3 - 14.4) / 17.3,NA)))))))))))))))

exp.data$gain.dif.tot <- ifelse(exp.data$party_combined==1 & exp.data$Treatment_polls==0, (6.8 - 8),
                            ifelse(exp.data$party_combined==2 & exp.data$Treatment_polls==0, (29.5-34.3),       
                                   ifelse(exp.data$party_combined==3 & exp.data$Treatment_polls==0, (6.1-7.4),
                                          ifelse(exp.data$party_combined==4 & exp.data$Treatment_polls==0, (5.4-6.3),
                                                 ifelse(exp.data$party_combined==5 & exp.data$Treatment_polls==0, (24.7-28.3),
                                                        ifelse(exp.data$party_combined==6 & exp.data$Treatment_polls==0, (3.1 - 4.3),
                                                               ifelse(exp.data$party_combined==7 & exp.data$Treatment_polls==0, (4.7-6.3),
                                                                      ifelse(exp.data$party_combined==8 & exp.data$Treatment_polls==0, (17.3 - 20),
                                                                             ifelse(exp.data$party_combined==1 & exp.data$Treatment_polls==1, (6.8 - 5.6),
                                                                                    ifelse(exp.data$party_combined==2 & exp.data$Treatment_polls==1, (29.5-24.8),       
                                                                                           ifelse(exp.data$party_combined==3 & exp.data$Treatment_polls==1, (6.1-4.6),
                                                                                                  ifelse(exp.data$party_combined==4 & exp.data$Treatment_polls==1, (5.4-4.8),
                                                                                                         ifelse(exp.data$party_combined==5 & exp.data$Treatment_polls==1, (24.7-20.9),
                                                                                                                ifelse(exp.data$party_combined==7 & exp.data$Treatment_polls==1, (4.7-3.5),
                                                                                                                       ifelse(exp.data$party_combined==8 & exp.data$Treatment_polls==1, (17.3 - 14.4),NA)))))))))))))))


exp.data$gain.dif.sq <- ifelse(exp.data$gain.dif<0,exp.data$gain.dif,exp.data$gain.dif^2) 
exp.data$gain.dif.pos <- ifelse(exp.data$gain.dif<0,exp.data$gain.dif*-1,exp.data$gain.dif)

exp.data$n.months.sq <- exp.data$n.months^2
exp.data$age <- 2016-exp.data$birth

exp.data$Next_election2 <- recode(exp.data$Next_election, "1=0;
                                                          2=0;
                                                          3=1;
                                                          4=2;
                                                          5=2")

# make variables predicting radicalization / moderation
rightwing <- c(1,4,5,6,8)
leftwing <- c(1,2,7)
exp.data$economy3 <- ifelse(exp.data$change_economy2==2 & exp.data$party_combined%in%rightwing | exp.data$change_economy2==0 & exp.data$party_combined%in%leftwing,2,
                            ifelse(exp.data$change_economy2==1,1,0))
# Omit non-treated parties
exp.data <- exp.data[-which(exp.data$party_combined>8),]
relevant.vars <- c("Satisfied_polls", "policy_change", "Change_strategy2", "economy3", "Next_election2","n.months", "sex", "age", "party_close_rec")
desc.data <- exp.data[,relevant.vars]
desc.data <- as.data.frame(stack(desc.data))

exp.data$policy_change <- factor(exp.data$policy_change)
exp.data$change_econ_rec <- factor(exp.data$change_econ_rec)
exp.data$change_welfare_rec <- factor(exp.data$change_welfare_rec)
exp.data$change_economy2 <- factor(exp.data$change_economy2)
exp.data$change_welfare2 <- factor(exp.data$change_welfare2)
exp.data$Change_strategy2 <- factor(exp.data$Change_strategy2)
exp.data$Next_election2 <- factor(exp.data$Next_election2)

exp.data$economy3 <- factor(exp.data$economy3)

# Do table 1 ----------------------------------------------------------
chisq.test(table(exp.data$Treatment_polls,exp.data$party)[,-2]) 
round((table(exp.data$party)/dim(exp.data)[1]) * 100,2)

var.names <- list("Satisfied_polls"="Satisfaction with polls",
"policy_change"="Change policy",
"Change_strategy2"="Change strategy",
"economy3"="Change position",
"Next_election2"="Prospect next election",
"n.months"="Number of months",
"sex"="Gender",
"age"="Age",
"party_close_rec"="Feel close to party")

labeller <- function(variable,value){
  return(var.names[value])
}

descriptive.hists <- ggplot(desc.data, aes(x=values)) +
  geom_bar() +
  facet_wrap(~ind, scales="free", ncol=3, labeller=labeller) +
  theme_minimal()

ggsave(descriptive.hists, file="figC1.jpg", dpi=900, width=5, height=5)  

# Do analysis treatment effect --------------------------------------------
t.test(exp.data$Satisfied_polls~exp.data$Treatment_polls)

model <- lm(Satisfied_polls ~ n.months + seem.to.lose + Treatment_polls, data=exp.data)
temp <- data.frame(mean(exp.data$n.months, na.rm=TRUE),expand.grid(c(0,1),c(0,1)))
colnames(temp) <- c("n.months", "seem.to.lose","Treatment_polls")
prediction <- as.data.frame(cbind(predict(model, temp, interval = "confidence"),temp))
prediction <- prediction[order(prediction$fit),]

prediction$label2 <- factor(c("Gains", "Losses", "Gains", "Losses"))
prediction$label <- factor(c("Recent gains", "Recent gains","Recent losses","Recent losses"))
#prediction$procedure <- factor(c(rep("Regression prediction\nwith covariates", 4), rep("Mean",2)), levels=c("Mean", "Regression prediction\nwith covariates"))

fig1 <- ggplot(prediction, aes(x=label2, y=fit, fill=label2)) +
  geom_bar(stat="identity") + 
  #ylim(0,5) +
  geom_errorbar(aes(ymin=lwr, ymax=upr), width=.5) +
  facet_wrap(~label, scales="free_x") +
  ylab("Level of dissatisfaction") + 
  scale_y_continuous(breaks=c(0,1,2,3,4,5), labels=c("", "Very satisfied", "Satisfied", "Neither", "Dissatisfied", "Very disatisfied"), limits=c(0,5)) +
  xlab("Treatment") +
  theme_minimal() + scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9")) + 
  theme(legend.position="none", 
        strip.background=element_rect(colour="#56B4E9", fill="#56B4E9"),
        strip.text.x=element_text(colour="white", face='bold'))

ggsave(fig1, file="fig1.jpg", dpi=900, width=5, height=5)  

covariates <- c("Months between poll and reference point", "Losing in recent polls", "Gains Treatment")
stargazer(model, type="html", dep.var.labels="Satisfaction with polls", title="Table C1. Regression coefficients predicting satisfication with polls", covariate.labels=covariates, out="table C1.htm")

# Do analysis policy change -----------------------------------------------

# Regression model with covariates
model1 <- clm(policy_change~ seem.to.lose + Treatment_polls + sex + age + party_close_rec, data=exp.data)
model1a <- lm(as.numeric(as.character(policy_change)) ~ seem.to.lose + Treatment_polls + sex + age + party_close_rec, data=exp.data)

model2 <- clm(change_econ_rec ~ seem.to.lose + Treatment_polls + sex + age + party_close_rec, data=exp.data)
model3 <- clm(change_welfare_rec ~ seem.to.lose + Treatment_polls + sex + age + party_close_rec, data=exp.data)

model4 <- multinom(change_economy2 ~ seem.to.lose + Treatment_polls + sex + age + party_close_rec, data=exp.data)
model4a <- multinom(economy3 ~ seem.to.lose + Treatment_polls + sex + age + party_close_rec, data=exp.data)
model5 <- multinom(change_welfare2 ~ seem.to.lose + Treatment_polls + sex + age + party_close_rec, data=exp.data) # does this effect disappear if you exclude soc.dem
model5a <- multinom(change_welfare2 ~ seem.to.lose + Treatment_polls + sex + age + party_close_rec, data=exp.data, subset=exp.data$party!="Social Democrats") 
model6 <- clm(Change_strategy2 ~ seem.to.lose + Treatment_polls + sex + age + party_close_rec, data=exp.data)

covariates <- c("Losing in recent polls", "Gains Treatment", "Gender", "Age", "Closeness to party")
stargazer(model1,model6, type="html", dep.var.labels=c("Change policy", "Change strategy"), title="Table C2. Regression coefficients predicting responses to polls", covariate.labels=covariates, out="table C2.htm")
stargazer(model4a,model5, type="html", dep.var.labels=c("Shift position", "", "Prospects next election", ""), column.labels=c("Stick","Radicalize","Same","Worse"), title="Table C3. Regression coefficients predicting responses to polls", covariate.labels=covariates, out="table C3.htm")


# Report models 1, 4, 5 & 6 
temp <- data.frame(expand.grid(c(0,1),c(0,1)),0,mean(exp.data$age, na.rm=TRUE),0)
colnames(temp) <- c( "seem.to.lose","Treatment_polls", "sex", "age", "party_close_rec")

predict.policy <- as.data.frame(predict(model1, temp, interval=TRUE))
predict.policy1a <- as.data.frame(predict(model1a, temp, interval="confidence"))

predict.policy <- predict.policy[,c(1,4,7)]
#predict.strategy <- as.data.frame(predict(model6, temp, interval=TRUE))
#predict.strategy <- predict.strategy[,c(1,4,7)]
#predictions <- rbind(predict.policy, predict.strategy)
#predictions$dep.vars <- rep(c("Pr(Stick to policy)","Pr(Stick to strategy)"),each=4)
predict.policy$losses <- rep(c("Recent gains","Recent losses"),2)
predict.policy$treatment <- c("Losses", "Losses", "Gains", "Gains")

fig2 <- ggplot(predict.policy, aes(x=treatment, y=fit.0, fill=as.factor(treatment))) +
  geom_bar(stat="identity") + 
  geom_errorbar(aes(ymin=lwr.0, ymax=upr.0), width=.5) +
  facet_grid(~losses) +
  #ylab("Level of dissatisfaction") + 
  #scale_x_continuous(breaks=c(0,1), labels=c("Gains", "Losses")) +
  xlab("") +
  ylab("Pr(Choosing policy status quo") +
  #ylim(c(0,0.5)) +
  theme_minimal() + scale_fill_manual(values=c("#999999", "#E69F00", "#999999", "#E69F00")) + 
  theme(legend.position="none", 
        strip.background=element_rect(colour="#56B4E9", fill="#56B4E9"),
        strip.text=element_text(colour="white", face='bold'),
        axis.text.x=element_text(size=10))

ggsave(fig2, file="fig2.jpg", dpi=900, width=5, height=5)  

predict.econ <- as.data.frame(predicts(model4a, "0-1;0;0;60;0"))
predict.econ$level <- factor(ifelse(predict.econ$level==0, "Moderate",
                      ifelse(predict.econ$level==1, "Stay", "Radicalize")),levels=c("Moderate","Stay","Radicalize"))
predict.econ$type <- rep(c("Recent gains", "Recent losses"),each=3)

fig3 <- ggplot(predict.econ, aes(x=level, y=mean, fill=level)) +
  geom_bar(stat="identity") +
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.5) +
  facet_wrap(~type) +
  #ylab("Level of dissatisfaction") + 
  xlab("Preferred policy position economy") +
  ylab("Probability") +
  ylim(c(0,1)) +
  theme_minimal() + scale_fill_manual(values=c("#999999", "#E69F00", "#999999", "#E69F00")) + 
  theme(legend.position="none", 
        strip.background=element_rect(colour="#56B4E9", fill="#56B4E9"),
        strip.text.x=element_text(colour="white", face='bold'),
        axis.text.x=element_text(size=10))

ggsave(fig3, file="fig3.jpg", dpi=900, width=5, height=5)  

direction.per.party <- ddply(exp.data, .(party), summarize,
                             mean=mean(Change_economy,na.rm=TRUE),
                             sd=sd(Change_economy, na.rm=TRUE),
                             se=sd / sqrt(length(Change_economy)))
### Hier was ik gebleven

polls <- read.csv("polls_sweden.csv",sep=";")
polls$month.year <- str_sub(as.character(polls[,1]),-8,-1)
polls$month.year[7] <- "Jul 2016" 
polls$month.year[11] <- "Jun 2016" 
polls$month.year[165] <- "Aug 2014"

polls.by.month <- stack(polls[,3:10])
polls.by.month$month.year <- rep(polls$month.year,8)
polls.by.month <- ddply(polls.by.month, .(ind, month.year), summarize,
                        mean=mean(values))  
polls.by.month$month.year <- as.yearmon(polls.by.month$month.year)
polls.by.month <- polls.by.month[order(polls.by.month$month.year),]
polls.by.month$y <- rep(seq(0,24,1),each=8)
polls.by.month$size <- ifelse(polls.by.month$mean<10,"Small parties","Large parties")
polls.by.month$ind <- recode(polls.by.month$ind, "'soc' = 'Social Democrats';
                                                  'mod' = 'Moderates';
                                                  'sd' = 'Sweden Democrats';
                                                  'mp' = 'Green Party';
                                                  'c' = 'Centre Party';
                                                  'v' = 'The Left';
                                                  'l' = 'Liberals';
                                                  'kd' = 'Christian Democrats'")
polls <- ggplot(polls.by.month, aes(x=y, y=mean, group=ind, colour=ind)) +
    geom_smooth() +
    facet_wrap(~size, scales="free") +
    scale_x_continuous(breaks=c(0,12,24),labels=c("Elections '14", "August '15", "August '16")) +
    ylab("Percentage votes in polls") +
    xlab("") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle=45)) +
    theme(legend.title=element_blank())

ggsave(polls, file="polls.jpg",dpi=900)
  
polls.by.month$position <- ifelse(polls.by.month$ind=="Social Democrats",3.89,
                           ifelse(polls.by.month$ind=="Moderates",7.94,
                           ifelse(polls.by.month$ind=="Sweden Democrats",8,
                           ifelse(polls.by.month$ind=="Green Party",3.12,
                           ifelse(polls.by.month$ind=="Centre Party",6.89,
                           ifelse(polls.by.month$ind=="The Left",1.56,
                           ifelse(polls.by.month$ind=="Liberals",6.56,
                           ifelse(polls.by.month$ind=="Christian Democrats",7.83,NA))))))))                                            

median.voter <- ddply(polls.by.month, .(y), summarize,
                      mv = sum(position*(mean/100)))

median.voter.plot <- ggplot(median.voter, aes(x=y, y=mv)) +
                     geom_line() +
                     scale_x_continuous(breaks=c(0,12,24),labels=c("Elections '14", "August '15", "August '16")) +
                     ylab("Centre of gravity") +
                     xlab("") +
                     ylim(3.5,6.5) +
                     theme_minimal() +
                     theme(axis.text.x = element_text(angle=45)) +
                     theme(legend.title=element_blank())
  

# Next election prediction ------------------------------------------------
model8 <- multinom(factor(Next_election2) ~ seem.to.lose + Treatment_polls + sex + age + party_close_rec, data=exp.data)

predict.elec <- as.data.frame(predicts(model8, "0-1;0;0;60;0"))
predict.elec$level <- factor(ifelse(predict.elec$level==0, "Better",
                                    ifelse(predict.elec$level==1, "Same", "Worse")),levels=c("Better","Same","Worse"))
predict.elec$type <- rep(c("Recent gains", "Recent losses"),each=3)

fig4 <- ggplot(predict.elec, aes(x=level, y=mean, fill=level)) +
  geom_bar(stat="identity") + 
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.5) +
  facet_wrap(~type) +
  #ylab("Level of dissatisfaction") + 
  xlab("Expectation next election") +
  ylab("Probability") +
  ylim(c(0,1)) +
  theme_minimal() + scale_fill_manual(values=c("#999999", "#E69F00", "#999999", "#E69F00")) + 
  theme(legend.position="none", 
        strip.background=element_rect(colour="#56B4E9", fill="#56B4E9"),
        strip.text.x=element_text(colour="white", face='bold'),
        axis.text.x=element_text(size=10))

ggsave(fig4, file="fig4.jpg", dpi=900, width=5, height=5)  


# Effects gov -------------------------------------------------------------
model1 <- polr(factor(policy_change)~ gov + seem.to.lose + Treatment_polls + sex + age + party_close_rec, data=exp.data)
model6 <- polr(factor(Change_strategy2) ~ gov + seem.to.lose + Treatment_polls + sex + age + party_close_rec, data=exp.data)


